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We present a general molecular-dynamics simulation scheme, based on the Nose thermostat, for 
sampling according to arbitrary phase space distributions. We formulate numerical methods based 
on both Nose-Hoover and Nose-Poincare thermostats for two specific classes of distributions; namely, 
those that are functions of the system Hamiltonian and those for which position and momentum are 
statistically independent. As an example, we propose a generalized variable temperature distribution 
that is designed to accelerate sampling in molecular systems. 



INTRODUCTION 

Molecular-dynamics (MD) computer simulation is a 
widely used tool in biology, chemistry, physics and ma- 
terials science 0,0. Much of the power in the technique 
lies in the ability to generate phase-space trajectories 
weighted according to a relevant statistical-mechanical 
distribution. In the first MD simulations, straightfor- 
ward integration of the equations of motion for the 
system under study yielded energy conserving trajecto- 
ries that, assuming ergodicity, generated microcanoni- 
cal (constant NVE) equilibrium distributions of phase- 
space configurations. Later, to better mimic exper- 
imental conditions, a variety of MD techniques were 
developed that generate phase-space distributions ac- 
cording to other standard statistical-mechanical distri- 
butions, such as canonical (NVT) || [|, || , isothermal- 
isobaric (7VPT)§, §, and grand-canonical(^VT)||. 
Recently, however, there has been growing interest in 
the simulation of systems with distributions that go 
beyond textbook statistical mechanical ensembles. For 
example, molecular-dynamics methods for the simula- 
tion of systems obeying Tsallis statistics] E| have been 
developed by Plastino and AnteneodoflC | and Andri- 
ciaoaei and Straub|ll[. In this work we outline a gen- 
eral molecular-dynamics scheme, based on the Nose 
thermostat H [t^, to generate configurations according 
to an arbitrary phase-space distribution. 

A primary motivation for the development of al- 
gorithms for the generation of non-standard distribu- 
tions is the need for methods that accelerate the con- 
figurational sampling of systems. Many systems are 
not sufficiently ergodic on the time scale of standard 
molecular-dynamics simulations to ensure the conver- 
gence of statistical averages. This is especially true of 
macromolecules, biomolecules and amorphous materi- 
als. Over the past decade, a number of methods have 



been developed to enhance sampling in MD. Berne 
and Straub|13 have recently written an excellent re- 
view of new sampling methods. Central to these ap- 
proaches has been the recognition that high activation 
barriers cause a bottleneck in phase space, rendering 
transitions between states unlikely. A common thread 
among many methods is the systematic deformation 
of the potential (or total) energy surface to accelerate 
barrier crossing, either by lowering the barriers or rais- 
ing the potential valleys. From a statistical mechanical 
perspective, such energy modifications induce a corre- 
sponding modification in the phase-space distribution 
by enhancing the statistical weight of configurations in 
the vicinity of energy barriers. Explicit knowledge of 
the modified sampling distribution allows for statisti- 
cal reweighting of the computed trajectories to achieve 
averages in the original ensemble. 

The simplest method for enhancing sampling is to 
scale the full Hamiltonian by some factor less than 
unity. This is equivalent to performing the simula- 
tion at a higher temperature. If averages are de- 
sired at temperature T, isothermal MD simulations 
can be carried out at some higher temperature T*, 
with averages at the original temperature computed by 
reweighting the probability of each configuration by a 
factor of exp [(-gj; — -gj^) H(p, <?)] ■ A significant dis- 
advantage with such temperature boost approaches is 
that, unless the boost is sufficiently small, low energy 
configurations are not visited with a frequency large 
enough to obtain acceptable statistics. A related ap- 
proach, Multicanonical MD, is based on a Monte Carlo 
technique of the same name) 14 and uses preliminary 
high temperature trajectories to construct a distribu- 
tion that is nearly flat in the position coordinates, al- 
lowing nearly uniform sampling of coordinate space in 
subsequent simulations. Multicanonical MD has been 
demonstrated to accelerate conformational sampling in 



model polypeptides |L5j and atomic clusters. 

Another approach is Voter's hyperdynamics Jl6|, |l7| ], 
which employs a "boost potential" to reduce the sam- 
pling probability in low energy regions, thereby accel- 
erating barrier crossing due to diminished relative en- 
ergetic cost. With boost potentials chosen to leave the 
potential energy in barrier regions unchanged, transi- 
tion state theory arguments can be used to obtain good 
transition rate statistics. However, the low-energy 
sampling problem remains as in high temperature dy- 
namics. Hyperdynamics has been used successfully in 
solid state systems, but the method is not generally ap- 
plicable to liquids, where the presence of many saddle 
points hampers the identification of well-defined bar- 
rier regions. 

Among recent methods of enhanced sampling, the 
approach of Straub and co-workers is of particular 
interest, with Monte Carlo and MD methods based on 
potential energy modifications that sample coordinates 
from alternative densities according to a formalism mo- 
tivated by the non-extensive Tsallis entropy ||. The 
Tsallis-Straub approach is easy to implement, amount- 
ing to a simple modification of the interaction forces 
according to the gradient of an effective potential. Re- 
cently, a more direct application of Tsallis entropy 
to MD was suggested by Plastino and Anteneodo|10|. 
Based on the idea of an effective Hamiltonian, these au- 
thors showed that canonical sampling with respect to 
the effective Hamiltonian is equivalent to Tsallis sam- 
pling. Significantly, this work considered only the Tsal- 
lis regime in which coordinate sampling is restricted to 
low energy regions. 

In this paper we present a dynamical framework 
for sampling according to a general class of proba- 
bility density functions, including but not limited to 
the Tsallis density. In order to introduce the idea of 
sampling from non-microcanonical distributions and 
to provide the necessary background for our gener- 
alized dynamics we discuss in Section the extended 
Hamiltonian approach of Nose|| to canonical (con- 
stant temperature) sampling, as well as the Nose- 
Hoover Q| and Nose-Poincare[^ approaches for imple- 
menting real-time formulations of Nose dynamics. In 
Section we introduce Generalized Distribution Dy- 
namics (GDD) and discuss the technique for two spe- 
cial classes of systems: those for which the position 
and momentum distributions are separable and those 
for which the phase space distribution is a function of 
the full Hamiltonian, and show how the Nose frame- 
work can be used to derive the equations of motion 
that produce trajectories that sample from generalized 
distributions. In section we present as an example 
the variable temperature distribution for accelerating 
the sampling of systems with high barriers. Numerical 
experiments on a double well potential using both the 
separable and full Hamiltonian GDD approaches are 
presented in section . 



SAMPLING FROM A CANONICAL 
DISTRIBUTION: THE NOSE THERMOSTAT 

In traditional (NVE) MD simulation, the equations 
of motion corresponding to the system Hamiltonian, 
i?(p,q), are integrated to generate the trajectories. 
The trajectory is constrained to the constant energy 
surface, E = H(p, q) determined by the initial values 
of coordinates and momenta. States in phase space 
along solutions are said to be sampled from the micro- 
canonical, or constant energy, distribution according to 
the probability density /?NVE(q, p) that is proportional 
to 5(H(q, p) — E), where S is the Dirac delta function. 

Due in part to a desire to bring simulation into ac- 
cord with laboratory experiments that are typically 
conducted at some fixed temperature, methods have 
been developed for generating trajectories which sam- 
ple from the canonical, or constant temperature, en- 
semble according to the probability density pNVT(q, p), 
which is proportional to exp [— (q, p)] where = 
l/(fcsT), T being the temperature and k B the Boltz- 
mann constant. In contrast to the microcanonical case, 
canonical sampling allows states at all energies, though 
higher energy states have lower probabilities depending 
on the value of temperature T. 

Although other methods exist, the most widely used 
techniques for generating canonically distributed tra- 
jectories in MD simulation are based on the extended 
Hamiltonian of Noseli, Ol: 
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where s and tt are conjugate thermostat variables, Q 
is a fictional thermostat mass which determines the 
strength of thermal coupling to the system, g — Nf + 1 
(with Nf being the number of degrees of freedom in the 
system) and p is a virtual momentum related to the 
actual momentum of the system by p = sp||. The 
equations of motion generated by the Nose Hamilto- 
nian (Eq. pi) are 
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The Nose method regulates the temperature of the 
system through a dynamical time transformation given 
by 4r — s, where t is the Nose (virtual) time and t is 
real time. The remarkable property of Nose dynamics 
is that microcanonical sampling of the extended phase 
space {q, p, s,7r} yields canonical sampling in the re- 
duced phase space, {q, p}, provided that the system is 
ergodic. 
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For practical calculations of averages such as veloc- 
ity autocorrelation functions, it is convenient to work 
in a real time implementation of the Nose thermostat. 
The most commonly used real-time modification is due 
to Hoover pj. Hoover recognized that one can generate 
a set of real-time equations of motion by making the 
following transformations to the Nose equations of mo- 
tion 

1. change of variables: p = p/s. 

2. time transformation: dr/dt — s, 

3. change of variables: rj = Ins and £ = f]. 

The result is the following time-reversible system 
of equations, known as the Nose-Hoover (NH) 
equations [Q: 

q = M x p (6) 
p = - W(q) - £p (7) 

V = £ (8) 
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[p T M- x p - gk B T] , 



(9) 



where g — Nf, the number of degrees of freedom 
in the system. These equations of motion are non- 
Hamiltonian in form since the coordinate transforma- 
tions were not canonical; however, a conserved energy 
does exist given by 
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gk B T V . (10) 



(Although the variable rj has been decoupled from the 
system, it is helpful to include it in the calculations so 
that E can be monitored as an indicator of trajectory 
stability.) 

Recently, Bond, Leimkuhler and Laird|| have devel- 
oped an alternative real-time Nose thermostat scheme, 
the Nose-Poincare method, which is Hamiltonian in 
form, allowing for the use of symplectic integration 
schemes (which have been shown to give superior sta- 
bility in long time simulation|]l8|). This is accom- 
plished by performing a time transformation, not to 
the Nose equations of motion as with Nose-Hoover, 
but directly to the Hamiltonian using a Poincare time 
transformation, as follows: 



H NP = s(H, 
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where Hq is the initial value of iJjvose- It can be easily 

verified^ that the phase space trajectory generated by 

Hn p is identical to that generated by if/vose except 

for a time transformation 4r = s. The Nose-Poincare 

at 

equations of motion are 
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Note that, the exact solution to Nose-Poincare equa- 
tions of motion generates trajectories that are identical 
to those generated by the Nose-Hoover scheme, exactly 
solved. It is in the construction of approximate numer- 
ical methods that these two approaches differ. 

Although we favor the Nose-Poincare method in 
all cases, the Nose-Hoover formalism is more familiar 
within the simulation community. For this reason, in 
the current article, we present schemes based on both 
the Nose-Hoover and Nose-Poincare approaches. 

In certain systems, for example those with few par- 
ticles or strong harmonic components, the ergodicity 
assumption basic to the Nose approaches is not met. 
For these cases, the notion of Nose-Hoover chains has 
been developed [jl9 1, in which the Hamiltonian is further 
extended with additional thermostat variables that are 
coupled to each other. It has been demonstrated that 
NH chains, with properly chosen thermostat masses, 
can induce the needed ergodicity so that NH dynam- 
ics provides a means of sampling from the canonical 
distribution. We discuss NH chains further in section 
and in the Appendix. 



GENERALIZED DISTRIBUTION DYNAMICS 

In this section we present a dynamical scheme for 
sampling points in phase space according to general 
function F(p, q) which satisfies the properties of a 
probability density function in the phase space vari- 
ables {p, q}: 



F(p,q) = landF(p,q)>0. 



In analogy to the procedure used by Plastino and 
AnteneodopJ to develop an MD method to generate 
the canonicaTTsallis distribution, we relate the general 
density to the canonical density by way of an effective 
Hamiltonian H c ff as 



-/3fleff 



which yields 



ff eff = --lnF(p,q) 



(17) 



It is clear that canonical sampling with respect to the 
effective Hamiltonian is equivalent to sampling accord- 
ing to the generalized probability density F. To achieve 
canonical sampling with H c s we write the Nose Hamil- 
tonian for Generalized Distribution Dynamics: 



H Nose = -^ln^(p/s,q) 



2Q 



gk B Thxs. (18) 
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From the equations of motion generated from this 
Nose Hamiltonian and after applying the transforma- 
tions described in the previous section, we obtain the 
Nose-Hoover GDD equations of motion: 
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Similarly, the Nose-Poincare equations of motion for 
GDD are 
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These formulations disrupt the separability of vari- 
ables present in the original NH and NP equations of 
motion [Eqs. and (13-16]), respectively]. A 

time reversible discretization of the GDD equations 
would involve the solution of nonlinear equations in 
q and p at every step. Iterative solution would re- 
quire many evaluations of the potential energy and its 
gradient at each step, likely adding tremendously to 
the computational burden. We address this issue by 
considering two special classes of probability density u 
functions that maintain variable separability: 

Case 1: GDD for Separable Distribution Functions 

Consider separable probability distribution func- 
tions of the form 

F(p,q L )=A(p)B(q). 

We can relate the separable density to the canonical 
density by way of effective kinetic and potential ener- 
gies K e g and V e g as 



F(p,q) =e- f>K *e- f,v « l , 



leading to 



K cS (p) = -^\nA(p): y eff (q)=-iln£?(q). (27) 

Canonical sampling with respect to the effective Hamil- 
tonian H e g = K e s + V e s is equivalent to sampling ac- 
cording to the generalized probability density F. Fol- 
lowing the procedure outlined in the previous section, 



canonical sampling with H c g c & n be achieved using the 
Nose-Hoover GDD equations of 

motion, which for a separable distribution function 
are obtained as: 



q = Vpifeff(p) 

P = -V q V r cff (q) - £p 

v = a 

e- 1 



[p T V p X cff (p) - gk B T] 



(28) 
(29) 
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Generation of the Nose-Poincare equations of motion 
for this class of distributions follows similarly. 

Note that these equations have a simple relationship 
with the NH equations @-(@). Any existing imple- 
mentation of the NH (or NP) equations of motion can 
be easily modified for separable GDD by the replace- 
ment of M _1 p by V-Keff(p) in equations (g), (||), and 
V(q) by Vcff(q) in equation (0). 

The most important applications for GDD for sepa- 
rable distributions arc those in which only the coordi- 
nate distribution is altered through modification of the 
potential. Such potential-only modifications are at the 
heart of Voter dynamics jL7| and the Tsallis statistics 
based methods for accelerated sampling of Straub and 
Andricioaeijll]]. For such systems K e g is equal to its 
standard form ip T M _1 p and V e g is given by eq. (p7|). 
Implementation of GDD for such systems is straightfor- 
ward as any existing Nose-Hoover (or Nose-Poincare) 
code could be used without modification (other than 
the use of a modified input potential surface). 

Case 2: GDD for distributions that are functions of the 
Hamiltonian 

Here we consider distributions that are formal func- 
tions of the scalar Hamiltonian: F{H{p, q)). Defining 
the effective Hamiltonian as 

H cS = {-1/(3) lnF(tf (p, q)) = f(H(p, q)) 

with associated Nose Hamiltonian, the Nose-Hoover 
GDD equations of motion for this case are 



q = /'(ff(p,q))M- 1 p 
p = -/'(ff(p,q))W(q)-£p 
f, = { 
1 



(32) 
(33) 
(34) 



i = ^ [f'(H(p,q))p T M- 1 p-gk B T] . (35) 

We can arrive at the natural expression for the mo- 
menta by performing the time transformation dt/dr = 
l//'(ff(q,p)): 



dq/df = M _1 p 
dp/dr = -W(q) 

drj/df 



W(p,q)) 



/Wp,q)) 



(36) 
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(38) 



4 



d^/dr = - 



1 



W(p,q)) 



gk B T 



(39) 



These equations have a suggestive form. The influence 
of the modified distribution is manifested solely in the 
thermostat variables. Deviation of f'(H) from unity 
can be viewed as a time-dependent scaling of simula- 
tion temperature T along with an inverse scaling of the 
thermostat mass Q. 

It is possible to rewrite these equations yet again 
to achieve separation of the coordinates and momenta. 
Just as the quantity in (fill) is constant along solutions 
of the NH equations, the uDD equations ((3^) (|3^) con- 
serve the related quantity 



£ / = /(ff(q,p)) + ^+.gfc B T7 ? 



Qe 



(40) 



Making use of our assumption that F and (and also/) 
is monotonic and hence one-to-one, we can solve (M): 



H(q, p) = r 1 I El - ^ gk B T V 



Qe 



and define a new function of 77 and £ 



4>M = vr (r 1 (eI 



gk B Ti] 



(41) 



so that the GDD equations become 



dqjdf = M x p (42) 

d P /dt = -w(q)-M,0£p (43) 

d V /dr = 0(ry,m (44) 

di/dr = i [p T M- x p - ^(r;,e). 9 fc s T] . (45) 



Equations 
thermostat varia 



introduce coupling between the 
ibles, but leave the coordinates and 
momenta separated, allowing for efficient discretiza- 
tion schemes. A similar trick can also be used to 
simplify the numerical calculations in the symplectic 
Nose-Poincare method. We discuss numerical methods 
in the appendix. 



EXAMPLE APPLICATION: A VARIABLE 
TEMPERATURE DISTRIBUTION 

The methods of this paper are very general in the 
sense that dynamical simulations can be made to sam- 
ple any smooth, invertible density function F(p, q). In 
this section we propose a particular distribution func- 
tion both for the purpose of demonstrating the meth- 
ods of this paper and to outline a potentially useful 
method for accelerating sampling in systems with high 
barriers. 

As mentioned earlier, one way to enhance the sam- 
pling of systems that are not ergodic on the time scale 



of standard simulation due to high barriers is to carry 
out the simulations at high temperature. The original 
distribution can the be recovered by reweighting the 
trajectory to compensate for the change in the distri- 
bution. However, for most situations the low energy 
configurations that have large weight in the original 
distribution are not sampled with sufficient frequency 
at high temperatures to yield adequate statistics after 
reweighting. 

To address the low energy sampling problem with 
high temperatures, we propose a generalized distribu- 
tion that has the effect of raising temperature only in 
high energy regions while leaving the low energy dy- 
namics unaffected. To begin, we define the monotonic 
function / 7 (s), designed to smoothly switch between 
the identity function and a linear function of slope 7: 



Ms) 




if s < so 

if so < s < si (46) 
if s > si 



where So < S < Si control the size of the switching 
window and the shape of the switching function, and 
the polynomial coefficients of the switching function 
are given by 



a 
b 



(2^-(l- 7 ) So -(l + 7) Sl )/( So -si) 3 

(2(l-7)^ + (2 + 7)( So si + S ?)- 

3*(*o + *i))/(*o -si f 

(7 s o(so + s si - 2s\) - 

si(4sg + s sx - 65s + sf))/(s - s x ) 3 

(sg<5(s - 3si) + 

si(2si + 7 (si - s ))) /(s - si) 3 . 



A graph of the function with 7 = 0.2, so = 3, Si = 4, 
and 8 = (so+si)/2is shown in the left view of Figure ||. 

We now define a probability distribution as a func- 
tion the Hamiltonian H : 



F(p,q) = F 7 (ff(p,q)) = e 



- „-/ T (H)/feT 



(47) 



The same switching function / 7 can be used to generate 
a separable distribution with 



F(p,q) = e 



_ p -p T M- 1 p/2fcT -/ T (y(q))/feT 



(48) 



in which the momentum distribution remains canon- 
ical. At 7 = 1 (and S = si) both the Hamiltonian 
and potential versions of the variable temperature dis- 
tribution reduce to the canonical distribution at tem- 
perature T. For 7 < 1, the distributions give canonical 
sampling in low energy (total or potential, respectively) 
regions (E < Eq for some predetermined value Eq) of 
phase space at reference temperature T while sampling 
high energy regions (E > E\) at the higher tempera- 
ture T(E) = T/a(E) where a(E) = 7(1 --§-) + -§• 
Note that a(E) approaches 7 in the limit of large 
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E. This distribution modification has similarities with 
the Tsallis-based distributions used by Andricioaei and 
Straubp| and Plastino and Anteneodo[jlO| in that 
the effective temperature is a monotonically increasing 
function of energy; however, in these cases the temper- 
ature (and thus the dynamics) is altered at all energies 
whereas in our present case the dynamics at low ener- 
gies is unaltered. 



NUMERICAL EXPERIMENTS 

To test our methods, we consider the double well 
potential 



V(x) = e(x 4 



2x 2 



1) 



with minima at x — ±1 and barrier height e. 

We have performed a number of GDD simulations 
using the variable temperature distribution (47) in 
both its full Hamiltonian and potential forms. For all 
simulations, the reference temperature was kT = j^. 
For such low dimensional systems it is necessary to 
use Nose-Hoover chains to enhance the ergodicity of 
the dynamics — see the Appendix for discussion. In 
all simulations six thermostats were used. The switch- 
ing window parameters in ( [46| ) were taken as so = 3, 
Si = 4 and S = 3.5 (except when 7 = 1, then S — si). 
We report experiments with the variable temperature 
parameter 7=1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2. 
As 7 decreases from unity, stability considerations dic- 
tate smaller timesteps. The timestep was h — 0.001 for 
the modified potential energy calculations. For the full 
Hamiltonian approach, we use timestep h = 0.0001. 

The left view of Figure fi] shows the distribution 
of coordinates for the full Hamiltonian calculations 
at 7 = 0.2. It can be seen by the good agreement 
with the theoretical distribution that the coordinates 
along the trajectory are sampled according to F 1 . Also 
shown is the reweighted distribution which recovers the 
canonical distribution. Note that at this temperature 
(kT = e/10) standard NH chain dynamics fails to sam- 
ple effectively. It must be pointed out that because 
the GDD equations for the full Hamiltonian case were 
derived with a time transformation dt/dr = </>(t7,£) in 
(41), it is necessary to include </> as a weighting function 
when computing averages using trajectories produced 
by equations (p2|)-(^5|). The right view of Figure [l] 
shows the distribution of coordinates for the modified 
potential calculations at 7 = 0.2. As for the results for 
the full Hamiltonian, it can be seen that the canonical 
coordinate distribution is also recovered by reweight- 
ing. Note that the unweighted coordinate distributions 
from the full Hamiltonian and modified potential for- 
mulations are not identical. In particular, the trajec- 
tory from the full Hamiltonian method spends slightly 
more time in high energy configurations. 

In Figure ^ we show the distribution of total energy 
along the computed full Hamiltonian trajectory for _F 7 






coordinate 



FIG. 1: Top view: Full Hamiltonian simulation p 7 with 7 = 
0.2, Bottom view: Modified Potential simulation with 7 = 
0.2. The solid curve gives the coordinate distribution of the 
computed trajectories, the dashed line gives the reweighting 
of the computed sampling to the canonical distribution, and 
the dotted lines give the theoretical canonical distribution. 



with 7 = 0.2, which can be seen to closely approximate 
the theoretical energy distribution for F~ . Also shown 
is the function f(H) from equation (Eq), along with 
the computed values of / calculated as the natural log- 
arithm of the energy distribution from the trajectory. 

In Figure || we illustrate the success of the variable 
temperature density _F 7 in hastening barrier crossings 
for the double well system. The figure shows waiting 
time plotted versus the temperature boost factor I/7. 
It can be seen that both the full Hamiltonian and mod- 
ified potential approaches yield dramatic reductions in 
waiting time between barrier crossings. 



CONCLUSION 

In this paper we have presented a general dynami- 
cal formalism, which we call Generalized Distribution 
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energy 

FIG. 2: Top View: Total energy density for the Full Hamil- 
tonian distribution; Bottom View: The function f(H) re- 
constructed from trajectory energy sampling. The solid 
curve gives the total energy distribution of the computed 
trajectories, the dashed line gives the reweighting of the 
computed sampling to the canonical distribution, and the 
dotted lines give the theoretical canonical distribution. 




High Energy Slope Reduction Factor y 



FIG. 3: Sampling speedup measured as increased frequency 
of barrier crossings for the full Hamiltonian (stars) and 
modified potential (circles) methods using the variable tem- 
perature distribution. The high energy slope reduction fac- 
tor is 7 _1 . The speedup factor is the average waiting time 
between barrier crossings for each value of 7, normalized 
by the average waiting time at 7 = 1. 



the method we performed numerical experiments using 
a one-dimensional bistable oscillator and demonstrate 
that the Variable Temperature Distribution is very ef- 
fective for accelerated sampling of coordinates when 
used in the effective potential energy setting. We are 
currently applying this work to enhance the dynamical 
sampling of model polypeptides. 
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Distribution, which has application in accelerating con- 
figurational sampling in systems with high energy bar- 
riers. To illustrate the numerical scheme and evaluate 
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Appendix 

Numerical Methods for Nose-Hoover GDD for 
distributions that are functions of the Hamiltonian 

The NH equations (|)-(|) can be discretized in a 
number of ways which generalize the basic Verlet ap- 
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proach. One such method is |2 
At 



At 

gk B T) (A-2) 

?7n+i = ??n+At^ n+ i (A-3) 

q n +i = Qn + AtM^Pn+j (A-4) 
£„j. i + — fp T j_ i M _1 p„, i - 
gk B T) 



+i 



At . 

Pn+i = P„+i - -y (VF(q„+i) + 



(A-5) 



(A-6) 



Equations ( |A-l| ) and ( |A-2| ) can be solved by computing 
the vector 



At 



P = Pr, 



and rewriting (A-2 



At 



: vy( q „) 



p^^p 

1 + lf£n+± 



This scalar equation for can now be solved analyt- 
ically or with an iterative solver. Using the computed 
£ n +i, we can compute 



Pn+i 



1 + ^rin+i 



The discretization scheme ( |A-l[) - 
eralized to the GDD equations J42|)-- 



can be gen- 



Pn+i = 



Vn+1 
Cn+1 



Pn - Y (VU(q„) + 

0(r? n ,^„+i)C„+|P„+i) 



4>(ri n ,£ n+ i)gk B T 



q„ + AtM- x p n+ i 
r?„+i + ^^(?7„+i,C„+i)C+i ( A " n ) 



(A-7) 

(A-8) 

(A-9) 
(A-10) 



At 
2Q 



Pn+1 



At 



Pn+A 



(W(q nH 



<Ktyi+l>Cn+i)ff»+AP7i+- 



(A-12) 



(A-13) 



Two of the formulae a re im plicit. The p rocedure for 



£ n+ 1 is the same as for ( A-8 ), while ( A-ll ) may require 



the use of an iterative method, depending on the nature 
of the function (p. In the numerical results presented 
here, we use Newton-Raphson iteration with tolerance 
of 10~ 12 in double precision. 

The GDD scaling function <j>(r),£) introduced in 
equation ( pll| ) serves two interesting purposes. The re- 
sulting dynamical formalism retains the form of the 
NH equations, with coordinates and momenta coupled 
to a generalized thermostat (subject to a time trans- 
formation) via the thermostat variables rj and £. From 
a practical point of view the introduction of (f> allows 
for discretization schemes which are explicit in the co- 
ordinates and momenta, which is important for overall 
efficiency . Im p lemen tation of a timestepping scheme 



such as ( A-7 )-( A-13 ) requires repeated evaluation of 
the function <f>, requiring the inversion of the func- 
tion / in which defines the effective Hamiltonian 
of the generalized density. In general an analytic ex- 
pression will not be available for For the work 
described here, we have implemented the GDD scal- 
ing function using an algorithm which relies on the 
Newton-Raphson method for finding a zero of a scalar 
equation. For the variable temperature distribution 
based on @, we evaluate <%£) = l/f'(f- 1 (H)) by 

-gk B Tr/ from equation 



first evaluating AE = E 3 Q 
j4l|), then evaluating 



oe 

2 



AE 

root ofas 3 + bs 2j r 
cs + d- AE 
^AE-S) + H t 



if AE < H 

if H <AE< H Y 
if AE > Hi 



and 




if A/- 1 < H 
if i^o < / 1 <#i 

if r 1 > #1 



With a good initial approximation (such as AE or the 
average ~(AE+±-(AE-5)+Hi)) the Newton-Raphson 
method above converges to the desired root with two 
or three iterations in our experience, subject to a con- 
vergence tolerance of 10~ 12 in double precision calcu- 
lations. 



8 



A symplectic numerical method for Nose-Poincare 
GDD for distributions that are functions of the 
Hamiltonian 

Starting from the Nose extended Hamiltonian ap- 
plied to the effective Hamiltonian, 



H 



Nc 



-f(H(q,p/s)) 



7T" 

2Q 



gk B T\ns (A-14) 



Assuming / is one to one, we can invert / to obtain, 
along the energy surface H^ osc = E, the new Hamilto- 
nian 



H(q,p/s)-f- 



I , f 7T 



gkTlns) = 0. (A-15) 



The zero energy dynamics in this Hamiltonian corre- 
spond to i?Nose = Eq dynamics. We next introduce a 
time-transformation of Poincare type, H — > sH result- 
ing in 



H 



f 



NP = S H(q,p/s)- S f- 1 (El 
2 





7T 



2Q~ 9kTlnS) 



(A-16) 



It is natural to use a splitting method here, break- 
ing the Hamiltonian into two parts according to the 
obvious additive decomposition and solving each term 
successively using an appropriate symplectic numeri- 
cal method. Note that the splitting suggested here is 
different than that used recently by Nose 0] in his 
variation of the Nose-Poincare method, but the basic 
technique is similar. Integration of the term 



H X = sH(q,p/s) 



(A-17) 



can be easily performed using the standard (and sym- 
plectic) Verlet method; note that during this fraction 
of the propagation timestep, s will be constant. For- 
mally, the integration of 



Nose-Hoover chains 

For systems that are small or contain stiff oscillatory 
components, lack of ergodicity may render the Nose 
scheme ineffective. Chains of Nose- Hoover thermostats 
have been shown []19| to allow canonical sampling in 
these cases. For a chain of m + 1 thermostat variables 
the equations of motion are 



dq/dr 
dp/dr 
drjo/dr 



M _1 p 

- W(q) - £oP 



(A-21) 
(A-22) 
(A-23) 



<%o/dT = — [p T M" 1 p - gk B T] - (A-24) 
Wo 

drfr/dT = £i (A-25) 

dii/dr = ±- [QoCo - k B T] - 66 (A-26) 



drim-i/dr 
d^m-i/dr 
dri m /dT 
d^m/dr 



£m— 1 
1 



(A-27) 



(A-29) 



(A-30) 



Along solutions of the Nose-Hoover Chain (NHC) equa- 
tions the conserved quantity is 



E 



P T M'p , M srQi& 

i=0 

m 

gk B Ti] + ^2 kBTrji. 

i=l 



(A-31) 



Discretization schemes for NH chains are discussed in 
References p2] and ]l9fl . Extension of Nose-Poincare 
to incorporate chains is somewhat delicate; this is work 
in progress by two of the authors. 



Ho 



2Q 



— gkT Ins) 



(A-18) 



can be done analytically. However, this is relatively 
painful. A simpler approach is to use an implicit 
method, such as the implicit midpoint method. For 
a general Hamiltonian H{q,p), the midpoint method 
advances from step to step by solving 

q«+i = q« + AiV pj ff(q I1+ i ,p n+ i) (A-19) 
Pn+i =p„ + AtV g #(q n+ i,p„ + i) (A-20) 

where q n+ i/ 2 = (q„ + q«+i)/2 and p n+ i/ 2 is defined 
similarly. In the present case, this means solving a 
nonlinear system in R 2 at each timestep. 
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